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We present a novel cluster-expansion (CE) approach for the first-principles modeling of tem- 
perature and concentration dependent alloy properties. While the standard CE method includes 
temperature effects only via the configurational entropy in Monte Carlo simulations, our strategy 
also covers the first-principles free energies of lattice vibrations. To this end, the effective cluster in- 
teractions of the CE have been rendered genuinely temperature dependent, so that they can include 
the vibrational free energies of the input structures. As a model system we use the phase-separating 
alloy Fe-Cu with our focus on the Fe-rich side. There, the solubility is derived from Monte Carlo 
simulations, whose precision had to be increased by averaging multiple CEs. We show that includ- 
ing the vibrational free energy is absolutely vital for the correct first-principles prediction of Cu 
solubility in the bcc Fe matrix: The solubility tremendously increases and is now in quantitative 
agreement with experimental findings. 

PACS numbers: 71.15.Mb, 71.15.Nc,81.30.Mh,63.20.dk 



First-principles modeling of phase stabilities of alloys 
is of scientific and technological importance. A major 
progress forward was made by the Cluster Expansion 
(CE) which is based on an Ising-like concept [ll-ll|. The 
power of CE consists in modeling concentration depen- 
dent properties of coherent alloy phases based on first- 
principles input information. For a system the energy 
Ece{o') for a particular atomic configuration a with N 
atoms is expanded in terms of hierarchical atomic ar- 
rangements such as points, pairs, triangles, and higher 
order objects. Those arrangements are called figures f, 
and the selected figure set is denoted by F. The CE then 
reads 



(1) 



f6F 



in which the geometrically determined correlations nt(CT) 
and the symmetry degeneracy D( are known for the given 
underlying parental crystal lattice. The unknown effec- 
tive cluster interaction energies (ECIs) Jf, which are in- 
dependent of (J, have to be extracted from some suitable 
input information, such as a set of density functional the- 
ory (DFT) structures, which arc denoted by cr G input. 
For those ordered structures the DFT calculations pro- 
vide the ground state total energies i?o.DFT(o'). Fitting 
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the CE to the DFT results determines the unknown Jf. 
This is performed by a least-squares minimization of the 
residuals [Bj, which in a simplified formulation 0,01 reads 



\Ece{o-) - £^o,dft(ct) 

crGinput 



■ mm 



(2) 



The fit is validated by a (leave one out) cross validation 
score (CVS) @, which in turn drives a genetic algorithm 
(GA) in order to select the optimal figure set F for the 
given input [1, 0, @l • Additional input is provided until 
the CE is converged in a self-consistent way. If the CE in 
Eq.[T] converges reasonably fast and the fit in Eq.[2]is suf- 
ficiently accurate then DFT accuracy can be carried over 
to a configuration space much larger than the one de- 
fined by the DFT input. Finally, the combination of CE 
with Monte Carlo (MC) simulations allows a tempera- 
ture dependent treatment of phase stabilities and related 
goperties for a very large number of interacting atoms 

So far, temperature only entered via the configura- 
tional entropy modeled by the MC simulation; other 
temperature dependent contributions were left out, e.g., 
the important vibrational free energies. In the following 
paragraphs, the present study will include the contribu- 
tions from lattice vibrations and will demonstrate their 
strong influence on the phase stability. 

Formally, it is obvious that the CE becomes tempera- 
ture dependent when the ECIs become temperature de- 
pendent: Jf Ji{T)- This is the result when the ECIs 



2 



in Eq. [2] are fitted to temperature dependent input en- 
ergies. In the present case those are obtained by sum- 
ming the temperature dependent vibrational free energy 
Fvih.DFT{<^,T) to the ground state total energy, 

EuFTicr, T) ^ Ea,T>FT{l7) + ^Vib,DFT(CT, T) , (3) 

in which i?o,DFT(o') is the outcome of a standard DFT 
calculation strictly valid only at T = OK. The label 
"DFT" for -Fvib.DFTlCT, T) indicates that it can be de- 
rived by the same DFT approach and accuracy as used for 
the total energy (see below for details). Other tempera- 
ture dependent properties may be included by adding the 
corresponding temperature dependent terms, such as the 
magnetic ordering energy. However, such contributions 
are not included in the present study and — regarding the 
magnetic ordering — we assume perfect ferromagnetic or- 
dering in terms of spin polarization. In order to include 
these temperature dependent effects, the CE is rewritten 
as 

Ecf{(t) -> EcE{a,T)^N ^ Di,h{T)Iii{a) . (4) 

feF(T) 

Note, that the optimal set of figures has also become 
temperature dependent: F — )■ F(r). 

In the following, we will perform and discuss the tem- 
perature dependent form of the CE where the addition- 
ally included vibrational free energy is in general im- 
portant for the phase stability of alloys and compounds 
[l0l - ll2| . For this purpose the phase separating binary 
Fei_a:Cua: alloy system at the Fe-rich side of the phase 
diagram is considered [l3j . For such a system the applica- 
tion of CE needs particular care because no ground state 
line of ordered compounds exists, i.e. all formation ener- 
gies are positive. Furthermore, besides the technological 
interest of hardening steel by alloying Fe with Cu, a pre- 
vious study based on isolated single-atom and pairwise 
defects indicated that vibrational free energies are indeed 
influential on the solubility of Cu in an Fe matrix 
Including vibrational contributions to CE has been previ- 
ously discussed [ll] and applied in very few cases [l6l.[T7j. 
The actual procedure, how to include the vibrational free 
energy is not unique. In the present work an approach is 
presented which — in combination with a fast and accu- 
rate procedure for deriving the phonon spectra — can be 
used in a convenient way for doing a CE and subsequent 
MC calculations. 

The DFT calculations for the total energies were done 
by the Vienna ab initio simulation package (VASP) with 
the pseudopotential construction according to the pro- 
jector augmented wave method p^l - l20j . The exchange- 
correlation functional was treated within the general- 
ized gradient approximation as parametrized by Perdew, 
Burke and Ernzerhof [2l[. All calculations were done 
spin polarized assuming ferromagnetic ordering of the 
Fe-atoms. Very good convergency of total energies and 
forces with respect to energy cutoffs and fc-point integra- 
tion was ensured. Accurate forces were derived for cal- 
culating the phonon spectra and vibrational free energies 
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FIG. 1: (color online) Enthalpy of formation derived from 
£'Ce/dft(c", r). DFT input values (various symbols) and CE 
predictions (black crosses) are compared. For the phonon cal- 
culations of each structure the percentage of imaginary fre- 
quencies is indicated. The random mixing energies are shown 
for T = K (standard CE, black dashed curve) and for 
T = 1200 K (CE with vibrational free energy; blue dashed 
curve) . 



by a direct force-constant method within the harmonic 
approximation as implemented in our program package 
fPHON, which is based on PHON 

All the CE and DFT calculations were made for Fe-Cu 
alloys with a bcc parental lattice, since the main interest 
is in the Fe-rich part of the phase diagram below the 
ferrite to austenite transition. For pure Cu, also the fee 
ground state total energy was calculated as a reference. 
For the CE the UNiversal CLuster Expansion (UNCLE) 
program package Q was applied. 

Initially, a standard CE for a bcc parental lattice was 
made utilizing only the DFT total energies for T = 
K. The results in Fig. [1] reveal that no stable binary 
phase for any composition exists, as it is expresse d by 
the positive formation energies. As expected [l^ Il4| . 
the configurations with the lowest formation enthalpies 
(and the form of the ground-state line) correspond to 
phase separating atomic arrangements, which consist of 
slabs of pure Cu and Fe. In total, an input DFT set of 
51 configurations was taken into account resulting in a 
CVS of 3.7 meV/atom at T = K. The input set in- 
cludes the energetically favorable structures as well as 
configurationally excited states in order to get reliable 
MC rcsuhs, cf. Ref. In Fig. [H we let the CE pre- 



dict the formation enthalpies of all 631 configurations a 
with unit cells up to 8 atoms large. The random mix- 
ing energy shown in Fig. [1] for T = K (no vibrational 
free e nerg y included) agrees well with the result of Liu 
et al. |24l |. With increasing temperature (i.e., including 
^vib,DFT(f , T) in the CE) the random energy is lowered 
and its maximum shifts to higher Cu concentrations, as 
shown in Fig. [U for T ^ 1200 K. 
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FIG. 2: (color online) Cross section through the 50 x 50 x 50 
Monte Carlo simulation cell (Fe atoms: black, Cu: red). The 
initial setup of pure Cu and Fe blocks as shown on the left 
panel is brought into thermodynamical equilibrium for a fixed 
temperature (right part). The volume in the Fe block, in 
which the dissolved Cu atoms are counted, is indicated by two 
green borders, which are three layers away from the interface. 
This ensures that no Cu atom of the Cu slab is erroneously 
counted as dissolved. The right most panel demonstrates the 
concentration of dissolved Cu solubility per layer with and 
without Fvib.DPT(r). CE and MC calculations were made for 
the merged figure set using averaged ECIs (see text). 



Including now i^vib,DFT(o', T) for all 51 structures. 
Fig. [T] reveals that a considerable number of configura- 
tions have phonon spectra with imaginary frequencies, 
which indicates dynamical instability. Since all configu- 
rations are not thermodynamically stable anyway (they 
have positive formation enthalpies), this is not surprising. 
Anharmonic coupling of phonon modes might possibly 
stabilize some of the phonon modes (25i] , but such a task 
is forbiddingly expensive. Therefore, the usual assump- 
tion of neglecting non-vibrating modes in the vibrational 
free energy is made. 

According to Eq. SI different temperatures yield dif- 
ferent F(r). However, one finds that the temperature 
dependence of the solubility is not as smooth a function 
of the temperature as expected. This is a direct effect 
of the GA [1, 0] selecting the figure set F(r), and it can 
indeed be likened to that kind of arbitrariness which en- 
ters even at a single temperature: n different runs of the 
GA yield n different ¥i{T). All of them are equally ca- 
pable to map the input data onto the CE (Eq. 2]) but 
yield slightly different results in MC simulations. For 
the usual CE applications, this does not pose a problem: 
the precision needed for MC simulations with respect to 
concentration is not as strict as needed here for the Cu 
solubility in Fe (< 1 at.%), as we will see later on. In 
our case we need a strategy which allows us both to find 
the expected smooth behavior of the solubility and to 
increase the precision of the prediction. 

We provide the following solution: an averaging pro- 
cedure either of the results (i.e., the solubilities) or — 
more physically — of the CEs themselves. For each tem- 
perature, n = 10 different CEs were constructed, with 
corresponding temperature dependent figure sets Fi(r) 



TABLE I: Results of 10 temperature dependent CE + MC 
runs and of one CE -f MC with the merged figure set (see 
text). Ni is the number of figures in the merged figure set 
F(T) (see Eq.fS}. The last two columns show the Cu solubility 
as an average value Xs(^i{T)) of 10 separate MC runs and as 
derived from averaged ECIs a;s(Fi(r)) (see Eq. [5]). 
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no Fvib,DPT(T): 


1150 


137 


0.19 ± 0.04 


0.18 


with Fvib,DFT(r): 


850 


130 


0.08 ± 0.03 


0.06 




1000 


125 


0.46 ± 0.10 


0.43 




1150 


118 


1.58 ± 0.23 


1.58 



and energies EcE.i{o',T). For each CE i, a separate MC 
run was performed, where the simulation took place in a 
50 X 50 X 50 supcrcell, starting with the phase separated 
system by dividing the MC cell into blocks of pure Fe 
and Cu (see Fig. [2]). This setup of fixed reservoirs of Cu 
and Fe atoms allows for an exchange of atoms between 
the two slabs using the Metropolis algorithm. Having 
reached thermal equilibrium at a given temperature, the 
solubility — i.e., the equilibrium concentration Xs{¥i(T)) 
of dissolved Cu which depends slightly on the figure set 
¥i{T) used — is determined by counting the dissolved Cu 
atoms in bulk Fe as sketched in Fig. [2l For the different 
CEs, the solubility scatters around the averaged value 
Xs{¥i) = J2^^i^s{Pi)/n. Table U shows in the column 
a;s(F,;(r)) that the fiuctuations become sizable at ele- 
vated temperatures because a high precision of the CE 
is needed to determine the Cu solubility at rather dilute 
concentrations. Therefore small fiuctuations of the CE 
have a significant impact on the solubility. 

Instead of running one MC simulation for each of the 
n CEs, the averaging scheme can also be applied to the 
CE sums. We note that averaging the results — i.e., de- 
termining Xs(¥i{T)) — is indeed different from averaging 
the CEs. The n single CEs (all with their own Fi(T) and, 
consequently, their own ECIs) are averaged: 



1 " 

Ece{<J,T) = -Y Ece4'^,T) =: N V DiJi{T)ni{a 

fGF(T) 



(5) 

On the right-hand side, we introduced the merged figure 
set ¥{T) = Fi(r) U • ■ • U F„(r) with its corresponding 
temperature dependent averaged ECIs Jf{T). Obviously, 
F(T) will comprise a larger number of figures (more than 
100 in the present case, see Table |l| than any individual 
CE (about 40 in the present case). It should be noted 
that the value of the ECIs Jt(T) is not the result of the 
CE fitting procedure in Eq. [2]but of the described merg- 
ing after the fitting. 

Table U compares the Cu solubilities averaged over n = 
10 MC runs with the result Xs{¥{T)) of one MC run using 
the merged figure set F(T) and the averaged ECIs of 
Eq. [SJ While both values agree very well within the error 
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FIG. 3: (color online) Phase boundaries of Fe-rich Fei-ajCui 
alloys. First-principles results of 10 CEs and one MC with- 
out (triangles, red line) and with vibrational free energies, 
utilizing temperature dependent ECIs (diamonds, red line) 
and merged figure sets; Shown are results averaged over 10 
corresponding MCs (red circles, dotted line) including error 
bars, and results of a first-principles calculation with single- 
atom and pairwise Cu-defects (blue dashed line) [T^]- Semi- 
empirical CALPHAD data [l^l are indicated as a solid black 
line . 

bars of Xs{Wi{T)), it is clear that the two approaches do 
not yield absolutely the same results, as already pointed 
out. 

Figure [3] presents the phase boundaries at the Fe-rich 
side. By comparing the results without and with contri- 
butions from the vibrational free energies the very strik- 
ing difference is obvious: without Fvib,DFT(o', T) the sol- 
ubility is much too small compared to semi-empirical 
CALPHAD data Obviously, vibrational entropies 

are responsible for this effect. A comparison of the CE 
+ MC derived phase boundaries to the isolated defect 
model reveals a perfect agreement at lower temper- 
atures. But at higher temperatures larger defect clusters 
of Cu atoms enter the stage, as demonstrated by Fig. S) 
The CE + MC simulation at 1000 K finds most of the dis- 
solved Cu as single-atom and pair-wise defects mirroring 
the isolated defect model. Increasing the temperature to 
1200 K, CE -I- MC produces a substantial percentage of 
larger sized Cu-clusters thus demonstrating the concen- 
tration dependence of this approach and the deficit of the 
isolated defect model. 

Summarized, we have presented a combination of CE 



and temperature dependent properties in terms of vibra- 
tional free energies. With the averaged CEs (Eq. [S]) a 
single set of ECIs Jf (T) within a merged figure set F(r) 
has been derived by which one can further study, for ex- 
ample, the growth kinetics of precipitates. The presented 
concept for a temperature dependent CE is in princi- 
ple straightforward and also feasible, in particular if the 
strategy of the merged figure sets is utilized. Clearly, 
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FIG. 4: (color online) Distribution of Cu cluster sizes given 
as percentage of the total number of dissolved Cu atoms for 
the point defect model [l3| (blue bars) and the temperature 
dependent CE -I- MC calculation (red bars) with the merged 
figure set strategy (see text). 



there is still need for future improvement: in particu- 
lar, one should aim at reducing the number of figures 
in the merged figure set in order to reduce the compu- 
tational cost of MC simulations. In the case of Fe-rich 
Fei_i:Cui: we have demonstrated that the inclusion of 
vibrational free energies in the CE-I-MC simulations is 
absolutely vital: only then are realistic values obtained 
for the solubility of Cu in an bcc-Fe matrix, and only then 
do our results agree with experimental data. The main 
physics behind this surprisingly large solubility of Cu in 
Fe is effectively described by a concentration and tem- 
perature dependent and purely first-principles approach 
which also includes vibrational free energies. 
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search Program VICOM (Vienna Computational Ma- 
terials Laboratory, project no. F4110). Calculations 
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